Environmental variables drive plant species composition and distribution in the moist temperate forests of Northwestern Himalaya, Pakistan

By assessing plant species composition and distribution in biodiversity hotspots influenced by environmental gradients, we greatly advance our understanding of the local plant community and how environmental factors are affecting these communities. This is a proxy for determining how climate change influences plant communities in mountainous regions ("space-for-time" substitution). We evaluated plant species composition and distribution, and how and which environmental variables drive the plant communities in moist temperate zone of Manoor valley of Northwestern Himalaya, Pakistan. During four consecutive years (2015–2018), we sampled 30 sampling sites, measuring 21 environmental variables, and recording all plant species present in an altitudinal variable range of 1932–3168 m.a.s.l. We used different multivariate analyses to identify potential plant communities, and to evaluate the relative importance of each environmental variable in the species composition and distribution. Finally, we also evaluated diversity patterns, by comparing diversity indices and beta diversity processes. We found that (i) the moist temperate zone in this region can be divided in four different major plant communities; (ii) each plant community has a specific set of environmental drivers; (iii) there is a significant variation in plant species composition between communities, in which six species contributed most to the plant composition dissimilarity; (iv) there is a significant difference of the four diversity indices between communities; and (v) community structure is twice more influenced by the spatial turnover of species than by the species loss. Overall, we showed that altitudinal gradients offer an important range of different environmental variables, highlighting the existence of micro-climates that drive the structure and composition of plant species in each micro-region. Each plant community along the altitudinal gradient is influenced by a set of environmental variables, which lead to the presence of indicator species in each micro-region.


Introduction
Mountains are the most remarkable landforms on the earth, representing different vegetation zones based on environmental variations [1]. They offer habitat heterogeneity based on microenvironmental variation along the altitudinal gradient [2,3], where environmental variables (including direct and indirect effects of abiotic and biotic effects) are important factors in determining altitudinal zone boundaries [4,5]. It is well recognized that altitude is a dynamic gradient along which several environmental variables [6] and species diversity [7,8] change concomitantly. Within this focus, plant biodiversity is strongly influenced by different environmental variables [9], and certain species can survive on the brink of extinction in high mountains across the world [10][11][12].
Many mountains across the globe are important hotspots of biodiversity [13][14][15][16], with roughly half of all plant species flourishing in hotspot areas [17,18]. However, despite this high endemism of species greatly influenced by various environmental gradients, such as edaphic, climatic and physiographic variables [19], these areas suffer a major impact from climate change [20,21]. Plant species present in these gradients have great adaptive power [22,23]; however, the speed with which climate change is advancing might not be sufficiently achieved by plant species adaptation in these areas, leading to a strong impact on the biodiversity of these areas [24][25][26], ultimately leading to variations in the community structure [27].
In the face of the current climate change and considering the importance of phytosociological studies for the understanding of biodiversity and species distribution, the Himalayas represent a fundamental piece for these studies. This region is facing a marked increase in temperature [28], which is three times greater than the global average [29]. This unprecedented rise in temperature, and modification of various environmental variables as well, may lead to shifts in species composition [24][25][26] and variations in community structure [27]. Many researchers have explored the effect of altitude on species composition, diversity, and forests formation structure during the previous two decades, with around half of these studies indicating an inverse association between species richness and altitude. Rahbek [30], on the other hand, did a quantitative investigation of altitudinal gradients of species richness and discovered that among plants, hump-shaped patterns of species diversity with peaks at mid-elevation are the most typically recorded, followed by monotonically declining patterns. Kluge et al.
[31] found a hump-shaped diversity pattern for seed plants in the Eastern Himalayas, even though endemic species richness peaks at higher elevations due to increasingly isolated habitats and smaller surface area in mountainous ecosystems, which promotes speciation [32]. Although there has been a considerable increase in the number of phytosociological studies in these altitudinal regions considered hotspots of biodiversity [7,19,22,23,33], there is still limited knowledge of how and which environmental variables drive the distribution and composition of plant species along altitudinal gradients in specific hotspots regions, such as the Northwestern Himalaya.
By assessing the patterns of composition and distribution of plant species in these biodiversity hotspots influenced by environmental gradients, we greatly advance our understanding of the local plant community and how environmental factors are affecting these communities, which is a proxy for assessing how impacts of climate change can affect plant communities located in mountainous regions [34]. In this context, we evaluated plant species composition and how and which environmental variables drive the plant species distribution of moist temperate zone of the Northwestern Himalaya, Pakistan. Briefly, we assessed (i) the potential plant communities present in the moist temperate zone; (ii) which are the environmental variables that most determine plant community structure in the moist temperate zone; (iii) whether there is variation in plant species composition between plant communities and which are the species that most contributed for species composition dissimilarity; (iv) whether there is variation of diversity indices among communities; and finally (v) which is the beta diversity process that most influence plant community structure in the moist temperate zone.

Study area
The present study was carried out in the moist temperate zone of Manoor valley (District Mansehra, Khyber Pakhtunkhwa), which is a mountainous valley (34.68165 N to 34.83869 N latitude, and 73.57520 E to 73.73182 E longitude Fig 1) in the Northwestern Himalayan belt of Pakistan [35][36][37][38][39] along an altitudinal range of 1932-3168 m.a.s.l. Monsoon winds are the main source of precipitation as well as the primary force of controlling erosion, topography, climate and vegetation of the Northwestern Himalaya [1].

Ethics approval and consent to participate
This study was approved by the Board of Study (BoS), Committee, Department of Botany and Advanced Study Research Board (ASRB) of Hazara University, Mansehra 21300, KP, Pakistan.

Vegetation sampling and plant identification
In different growing seasons (from March to October), we evaluated 30 sampling sites per year, during four consecutive years, from 2015 to 2018. The line transect method (50 meters) we used for quantitative samplings [40][41][42][43][44][45], but we never repeated the same transects over years. The surveyed study area was subdivided into 30 stands and three points randomly selected were sampled within each sampling site along 50 meters transect (total = 90 transects).
The distance between the sampling sites was kept at 200 meters, while the distance between the transects was kept at 100 meters. The individuals of plant species falling precisely on the line were recorded. The data from each sampling site was calculated using phytosociological characters (i.e., density, frequency, cover and their relative values, as well as importance value) [46][47][48]. The IV was further used to rank each plant species and species with the highest IV were considered as the dominant species [46,47]. Similarly, each plant community was named based on three dominant species [49][50][51][52]. The slope angle, aspect and exposure were recorded using clinometer; and altitude, longitude and latitude were measured by geographical positioning system (GPS). Plants species collection, labelling, pressing and other herbarium work methodology was adopted following Ijaz [53], Ijaz et al. [54], Amjad et al. [55] and Stefanaki et al. [56]. Identification was done with the aid of Flora of Pakistan [57][58][59] and submitted to the Herbarium of Hazara University, Mansehra (Pakistan).

Environmental gradients
Soil samples weighting 200 grams were collected (15-30cm depth) randomly from each transect of sampled vegetation site [60,61]. The replicated samples of each sampling site were thoroughly mixed to form a composite sample [62], which was placed in a sterile polythene bag and labelled accordingly. Raw materials such as rocks, and stones were sieved out and the samples were then shade dried. Each dried soil sample was processed for physicochemical tests [62,63] to determine soil texture (i.e., clay, sand, silt, loam), pH [64], electrical conductivity (EC) [65], organic matter (OM) [66], nitrogen (N), phosphorus (P), potassium (K) and calcium carbonate (CaCO 3 ) levels [60,61,67]. Other climatic variables such as barometric pressure, dew point, humidity, heat index, temperature, wet bulb (relative humidity and ambient air temperature) and wind speed were also determined using a small remote weather station (Kestrel weather tracker 4000) to record the data at each transect and then average values were calculated at sampling site level [19].

Statistical analyses
The recorded data of vegetation, edaphic, and other environmental variables of sampling sites were compiled in order to determine relationship among them [68,69]. Matrixes of IV data of all the recorded plant species (244 species) from 30 sampling sites were used in the analyses. Analyses were conducted on PC-ORD [70] and R software 4.0.0 [71]. Packages used in R software are mentioned in each analysis. A georeferenced map was prepared to show the distribution pattern of distinct plant communities (Fig 1).
Sampling effort and cluster analyses. Species area curves (SAC) were used to check the efficiency of the sampling effort, where plant abundance data with Sørensen distance values were used to create species area curves [72]. For classification of the recorded plant species (244) and 30 sampling sites into different plant communities, we used three different cluster analyses: Two-way Indicator Species Analysis (TWINSPAN), and Cluster Analysis (CA). We identified and classified plant species and stands (sampling sites) into major plant communities [73], as well as assess the effects of various factors (such as environmental variables) on such communities by processing clustering method using TWINSPAN [19,74].
Plant communities and associated environmental variables. Both species and stands (sampling sites) were constrained in relation with the environmental gradients [75,76], which were divided into geographic, slope aspect, edaphic and climatic gradients. We used non-metric multidimensional scaling (NMDS) and Principal Component Analysis (PCA) ordination biplot to determine the relationship between vegetation communities and environmental gradients using the "vegan" package. In NMDS and PCA, arrows represent the environmental gradients, in which the length shows the strength of the gradient, and the direction represent the degree of correlation. The direction of gradients on the same axis reveals a positive correlation.
In addition, we performed canonical correspondence analysis (CCA) and variation partitioning tests (partial CCA) [77] to observe how explanatory variables (climatic, edaphic, geographic, and slope) drive the plant species distribution. First we built the best model with the lowest number of variables (those that most explain variance), through the step function with permutation using the "stats" package [71]. Next, we also evaluated multicollinearity between variables of the final model using Variance Inflation Factor (VIF), and we removed any variable with VIF>10, one at a time. Finally, with the final model, we carried out CCA and partial CCA to check how much each group of variables (geographic, edaphic, climatic, slope) explain in our model [77]. For these analysis, we used the "vegan" package [78].
Variation of plant species composition among plant communities. To compare whether there is difference in species composition between plant communities, we also used NMDS followed by a Permutational Multivariate Analysis of Variance (PERMANOVA) with Euclidean distance and 999 permutations. After PERMANOVA, we conducted pairwise comparisons between communities with corrections for multiple testing also using Euclidean distance and 999 permutations. We used false discovery rate (FDR) as p-value adjustment method. PERMANOVA and pairwise comparisons were conducted with "RVAideMemoire" package [79]. To observe the contribution of each plant species to overall dissimilarities, we used a similarity percentage analysis based on the decomposition of Bray-Curtis dissimilarity using the package "vegan" [78].
Variation of diversity indices among plant communities. To compare the diversity indices evaluated (species richness, Shannon, Simpson, and Pielou) among the four communities, we also conducted a GLM with Gaussian error distribution, except for species richness, in which we used Poisson distribution followed by Likelihood Ratio test. Pairwise comparisons were conducted with estimated marginal means using the package 'emmeans' [80].
Beta diversity. To evaluate which is the beta diversity process that most influence plant community structure in the moist temperate zone, we decomposed the Sørensen dissimilarity index (βsor), a measure of overall species replacement into two additive components: the spatial turnover (Simpson pairwise dissimilarity, βsim) and nestedness-resultant components (nestedness-fraction of Sorensen pairwise dissimilarity, βsne) [81][82][83]. Dissimilarity analysis was conducted in the package "betapart" [84].

Sampling effort and plant communities
In total, 244 plant species belonging to 194 genera and 74 families (Table 1) were recorded in the moist temperate forests of the Manoor valley, Northwestern Himalaya, Pakistan. The moist temperate forests ranged from 1932.3m to 3168m. SAC analysis showed that the maximum number of plant species appeared up to stand 26 and the species curve became parallel after it, as no new species were recorded further (Fig 2).

Plant communities and associated environmental variables
NMDS and PCA were used to show the relationship between the plant communities of moist temperate forests and environmental variables (Fig 4A-4D) and PCA (Fig 4E). The ecological and environmental variables like geographic, slope, edaphic, and climatic variables were used to correlate communtities ( Table 2). The most representative environmental variables that drive the community structure and diversity were altitude, slope angle and aspects (SE, NE, ES, WN), potassium (K), pH, organic matter, loam, silt, sand, clay, temperature, heat index, wind speed and barometric pressure. Environmental variables classify 30 sampling sites into four major plant communities, as shown by the cluster analysis (Fig 3). In constrained PCA ordination, the PC1 axis accounted for the most explanatory variance (20%), while the PC2 axis accounted for the least (14.2%). The profound influence of the environmental variables was revealed by classifying the moist temperate forests vegetation into four communities (Fig 4E), as also shown by CA, TWCA and NMDS.
The PCP community showed positive and significant correlation with northern aspect, silty loamy soil texture, humidity and altitude (Fig 4B and 4C). In contrast, CPI community showed positively significance with southern slope, wind speed, dew point and wet bulb. IHC community showed positive correlation with silty soil texture, electric conductivity and pH. And   Table 2. IHC: Indigofera heterantha-Heracleum candicans-Cynodon finally, VPI community revealed positively significant correlation with north-western slope, CaCO 3 and pH. Thus, all the four communities were found separately in clumps with clear differences based on the environmental variables (Fig 4).
The CCA and variation partitioning tests showed that the total inertia results of CCA was 3.023, where our final variables (altitude, temperature, humidity, wind speed, slope angle, slope N, slope NW, slope SW, pH, EC, OM, CaCO 3 , K, P, sand, and loam) together explained 66

Variation of plant species composition among plant communities and beta diversity
We found a significant variation in plant species composition among communities ( Table 5; Fig 6), in which all communities showed a significant difference in species composition between each other ( Table 6). Out of 244 species, six species greatly contributed to the variation in plant species composition between communities, namely Viburnum grandiflorum, Indigofera heterantha, Heracleum candicans, Cedrus deodara, Pinus wallichiana, and Parrotiopsis jacquemontiana (Table 6). Overall, the three species that most contributed for the variation in species composition between communities showed 13.7-29.7% of cumulative contribution ( Table 6). The total beta diversity (βsor) showed a value of 54.7% dissimilarity, of which spatial turnover (βsim) made up 40.5% and nestedness-resultant components (βsne) made up 14.2%. In βsim cluster, we observed 47.8% dissimilarity between PCP-CPI cluster and VIP-IHC cluster (Fig 7). PCP showed a dissimilarity of 9.4% with CPI, and VIP showed a dissimilarity of 21.5% with IHC (Fig 7). In βsne cluster, we found 24.5% dissimilarity between PCP and VIP-C-PI-IHC cluster (Fig 7). VIP showed a dissimilarity of 11.3% with IHC-CPI, and IHC had 4.1% Table 3. The contribution and ranking of the studied variables in the variation partitioning tests (partial CCA model) to observe how explanatory variables (i.e., climatic, edaphic, geographic, and slope) drive the plant species distribution. dissimilarity with CPI (Fig 7). Thus, plant community structure is twice more influenced by the spatial turnover of species (βsim) than by the species loss (nestedness-resultant, βsne).

Discussion
Mountain ecosystems are characterized by dramatic changes in temperature and abiotic properties over short altitudinal and geographical distances, making them ideal natural laboratories for studying vegetation response to environmental changes [85]. In this study, we evaluated the plant species composition and distribution in a hotspot of biodiversity, the Northwestern Himalayan mountains, Pakistan, assessing how environmental gradients, source of habitat heterogeneity, influence plant community structure and diversity, which might be a proxy for assessing how climate change impacts on plant communities located in mountainous regions Table 5. PERMANOVA results comparing species composition between the four communities found in Moist temperate forest. This analysis was made with Euclidean distance and 999 permutations. Pairwise comparisons between communities are depicted in Table 6.  Table 6. Pairwise comparisons with FDR p-value adjustment method of species composition and contrast results of the contribution of individual species to the overall Bray-Curtis dissimilarity of species composition between the four communities found in moist temperate forest. We displayed only the three species that most contributed. [34, 86,87]. We found that (i) the moist temperate zone in this region can be divided in four different major plant communities; (ii) each plant community has a specific set of environmental drivers; (iii) there is a significant variation in plant species composition between communities, in which six species contributed most to the plant composition dissimilarity; (iv) there is a significant difference of the four diversity indices (species richness, Shannon, Simpson, Pielou) between communities; and finally (v) plant community structure is twice more influenced by the spatial turnover of species (βsim) than by the species loss (nestedness-resultant, βsne). Overall, we showed that altitudinal gradients offer an important range of different environmental variables, highlighting the existence of micro-climates that drive the structure and composition of plant species in each micro-region. In addition, each plant community along the altitudinal gradient has a set of environmental drivers, which lead to the presence of indicator species in each micro-region. Mountain plant communities are thought to be sensitive to climate change and, thus, able to reveal its effects sooner than others [34,88]. The four communities found showed a wide range of environmental drivers; however, altitude and temperature showed great prominence, probably making up the main environmental drivers in mountainous plant communities. Similar pattern was observed in the allied area (Nandiar catchment, Battagram) of Northwestern Himalaya by stating altitude and temperature as the governing gradient [74]. Such variables, which can be strongly correlated [89], modify the diversity and structure of plant communities by creating local micro-climates [90], directly influencing plant community composition and diversity [19,26,91,92].

P-value
Indeed, there is no order of importance of environmental variables, but studies are unanimous in showing that there is a consensus on the explanations for the variables' influences. For instance, in two recent studies we showed that the altitude-temperature relationship significantly influenced the physiological attributes of some plant species in the Northwestern Himalayan region [22,23], which can be a proxy for understanding plant adaptation to climate change. Any change in soil parameters has a significant effect on the growth of plant communities [19]. The studies on mountain forests habitats around the world have also revealed the role of soil structure on species zonation [72,93,94]. Furthermore, both chemical and physical attributes of the soil are related to natural soil characteristics, with an impact on plant species composition and distribution of higher vascular plants [95][96][97]. For instance, some soil variables can have great influence on plant composition and distribution, such as pH. Some studies have shown that pH level on soil can influence nutrient availability, ultimately influencing nutrient uptake for growth [98][99][100]. However, the availability of some nutrients as a result of pH levels can be detrimental for some plants, since some nutrients are toxic to some plants [98,101]. Considering that there is a great variability of pH levels and nutrient availability and concentration along altitudinal variables, it is expected a great variability of plant species composition which can be more or less related to specific soil parameters. Since plants are sensitive to small variations of soil characteristics such as pH, minerals, organic matter, among others, and these variables are constantly changing along altitudinal gradients directly and indirectly influencing the presence and availability of other organisms and resources, some plant species might have adapted to specific set of variables.
Variability in plant species diversity is an outcome of species interaction with particular set of environment variables either abiotic and biotic [102,103], which can occur in both space and time [104,105]. The concept of changing species composition and vegetation continuum along the ecological gradients emerged as an antithesis model for distinct units [106,107]. In our study, the moist temperate forest of the studied Northwestern Himalayan region is comprised in an altitudinal gradient of approximately 1500 m. This gradient is subject to strong micro-climatic variation, which results in a set of micro-regions (better discussed above). Each micro-region has certain characteristics, which will influence the set of species that will inhabit these spaces [24][25][26]. In this sense, it is expected that the plant community structure is more influenced by the spatial turnover of species (βsim) than by the species loss (nestedness-resultant, βsne), i.e., that there are different plant communities along the altitudinal gradient, as shown by our results. The differentiation of species diversity was mainly a consequence of environmental variables which is due to soil factors [108]. Therefore, in addition to the influence of edaphic factors in space on species composition and vegetation continuum, as shown in our study, results from similar studies have shown that the altitude is also important in driving vegetation structure and diversity in plant communities.
We found that environmental heterogeneity among plant communities have significant effects on beta diversity, particularly the spatial turnover. These results indicate that there is not a significant loss of the number of species between the plant community, but a variation in the species composition. This variation may be closely linked to the environmental effects in the area, which induces the appearance of species adapted to environmental variables [109]. The local community composition replacement implied the simultaneous loss and gain of species due to immigration-extinction dynamics and trait-based environmental filtering [110,111]. This indicates the relationship among plant community types and among species based on multiple factors. Although we did not find a large variation in βsne (loss of species between plant communities), it is important to note that temporal analysis might be important to consider a notable variation in this component of beta diversity; and βsne variations will be better observed in long-term analysis in future studies. In this sense, it is expected that the plant community structure is more influenced by the spatial turnover of species (βsim) than by the species loss (nestedness-resultant, βsne), i.e., that there are different plant communities along the altitudinal gradient, as shown by our results. Similarly, results were report by Haq et al. [112] from forests of Kashmir Himalaya, India.
We observed a significant variation in plant species composition between communities, in which all communities showed a significant difference in species composition between each other. The measure of Bray-Curtis dissimilarity shows that species composition change that is influenced mainly by abundant species, in our study six species (Viburnum grandiflorum, Indigofera heterantha, Heracleum candicans, Cedrus deodara, Pinus wallichiana, and Parrotiopsis jacquemontiana) contributed most to the plant composition dissimilarity. These results suggest that the richness and turnover patterns we observed were driven primarily by rare species, which comprise most of the local species pools at these forest communities [113]. These findings are consistent with the idea that less abundant species are more sensitive to climate variability than longer lived and more abundant species [114]. The high level of turnover is common and is an important mechanism by which a large regional species pool buffers site level diversity from interannual variation in climate [115].
Current study provides the baseline and first insights of spatial distribution, vegetation pattern and species contribution in response to environmental gradients in a moist temperate forests, Northwestern Himalaya, Pakistan. Studies that evaluate the distribution and composition of the plant community are fundamental for a better understanding of the local plant community, the conservation status and protection of these communities, as well as providing support for mitigation measures. Especially in the case of Northwestern Himalaya, which represents a biodiversity hotspot, it is even more important that we conduct phytosociological studies in these areas to document and preserve the biodiversity there. In the face of current climate changes, these regions are being heavily impacted [28,29], where the probability of species extinction may be higher than elsewhere, as these regions are rich in endemic species. Finally, we need to consider that phytosociological studies consider a general profile of the first trophic chains level, i.e., to evaluate the composition, distribution and diversity of plants is to indirectly assess the first level of trophic chains.